\documentclass{article}
\usepackage{ctex}
\usepackage{amsmath}
\usepackage{appendix}
\usepackage{subfigure}
\usepackage{epstopdf}
\usepackage{titletoc}
\usepackage{titlesec}
\usepackage{pdfpages}
\usepackage[boxed,linesnumbered]{algorithm2e}
\usepackage{geometry}
\usepackage{graphicx}
\usepackage{listings}
\usepackage{float}
\usepackage{tikz}
\numberwithin{equation}{section}
\numberwithin{figure}{section}
\geometry{top=3.5cm,bottom=3cm,left=1cm,right=1cm}
\linespread{1.25}
\titleformat{\section}{\centering \zihao{4}}{第\thesection 章}{2em}{}
\titleformat{\subsection}{}{\thesubsection}{2em}{}
\renewcommand{\contentsname}{\centering 目录}
\newcounter{dirary}
\setcounter{dirary}{21}
\begin{document}
	\begin{titlepage}
		\includepdf{fpage.pdf}
	\end{titlepage}
	\dottedcontents{section}[2em]{\large}{2.0em}{12pt}
	\tableofcontents
	\thispagestyle{empty} %暂时取消页码
	\newpage
	\songti \zihao{-4} 
	\setcounter{page}{1}  %设置页码开始
	\section{实训目的和意义}
	%地球物理勘探是利用物理学的原理和方法，对地球的各种物理场的分布及其变化进行观测，并对观测的物理场进行处理和分析，获得地球内部及近地空间的构造、介质组分、物理空间的形成和演化历史等信息，并据此研究各种自然规律和矿产资源分布的一门学科。地球物理勘探中最核心的部分包括数据的采集、处理和解释。用地球物理方法研究或勘查地质体或地质构造，是根据测量数据或所观测的地球物理场求解场源体的问题，是地球物理场的反演的问题。为了研究地球物理场的传播、分布规律或验证处理方法是否正确，也需要对物理过程进行数值模拟。这些都离不开计算机编程，通过将推导的数学物理方程编写成计算机可以处理的程序，对波场进行正演或对数据进行处理。因此，作为地球物理类的学生，熟练掌握一门计算机语言并具备利用计算机编程解决实际问题的能力十分重要。本实训通过信号的生成、一维和二维傅里叶变换、声波方程数值模拟等内容的练习，使学生熟练掌握文件的读写、指针、高维数组及函数等的综合运用，提高他们利用计算机和数学物理方法解决实际问题的能力。
	作为一名物探本科生，学习的专业知识多为探究地球各种物理场。为了更加直观的比较地球各种物理场的差异，我们往往需要利用所测得的参数进行大量数据计算和分析。在数据分析的过程中，不可避免的要使用计算机来进行编程。C语言具有移植性好、执行效率高、占用内存小且简单易学等特点，因此使用C语言编程不仅解放了我们的双手，而且经济高效。
	
	本次实训是对我们之前所学C语言理论知识的一次实践应用，实训内容包括地震子波的生成、一维二维离散傅里叶正逆变换以及波动方程数值模拟，这些内容贴近地球物理实际数据处理过程，锻炼了我们动手解决问题的能力。此外，训练中多次用到对文件的读写操作和动态内存的开辟，训练了我们使用C语言进行文件操作和使用指针的能力，对于计算分析大量数据有很大帮助。
	
	本程序是在Linux平台上开发，因为不能使用Fimage软件，所以图像是由Seismic Unix绘制完成。
	\section{雷克子波}
	\subsection{准备知识}
	实际生活中，地震子波是一个很复杂的问题，因为地震子波与地层岩石性质有关，而地层岩石的性质本身就是一个复杂体。为了方便模拟地震子波，人们对地震子波进行了简化。目前普遍认为雷克提出的地震子波数学模型具有广泛的代表性，这个子波即为“雷克子波”。其表达式如下
	\begin{equation}
		\label{rick}
		f\left(t\right)=[1-2(\pi ft)^2]e^{-\left(\pi f t\right)^2}
	\end{equation}
	
	在上式中，$f$表示主频，$t$表示时刻。
		\subsection{程序设计}
	为了在时间的正半轴得到雷克子波的完整图像，我们可将子波沿时间轴向右平移一个周期的长度。根据式\ref{rick}，我们可将平移后的子波函数写成如下形式
	
	\begin{equation}
		\label{newrick}
		f\left(t\right)=\{1-2[\pi f(t-T)]^2\}e^{-\left[\pi f (t-T)\right]^2}
	\end{equation}
	
	其中$T$表示子波的周期，其值可表示为主频的倒数。
	
	式\ref{newrick}可以用下面的编程语句表示出来。其中，$fpeak$表示雷克子波的主频，即式\ref{newrick}中的$f$；$t$表示时刻，即式\ref{newrick}中的$t$；$T$代表子波的周期，即式\ref{newrick}中的$T$。
	\newpage
	\begin{algorithm}
		\SetAlgoNoLine
		\SetKw{pie}{define PI 3.1415926\\}
		\SetKwBlock{for}{for(it=0;it<nt;it++)\{}{\}}
		\pie
		T=1./fpeak\;
		\for{t=it*dt-T\;
			ric[it]=(1-2*PI*fpeak*t)*exp(-(PI*fpeak*t)*(PI*fpeak*t))\;}
	\end{algorithm}
	
	程序编写过程中，雷克子波取样点的数量是由用户输入而确定的，因此我们并不能在程序运行之前确定代表子波的数组的长度大小。我们只能在程序编写过程中为数组分配足够大的内存或者是采取动态分配的方法。若使用第一种方法，当用户输入的取样点数目过少，就会造成内存的浪费，而取样点过多以至于超过事先分配的内存，又会造成数据溢出。因此我们选择动态分配的方法。但动态开辟内存时，并不是每次都会如愿，若开辟失败则会使程序运行失败或访问到错误的内存地址，从而得到错误的数据，导致程序失败。为了确保程序正确运行，避免动态分配失败对程序的运算结果造成影响，我们需要知道动态分配成功与否。
	
	\begin{algorithm}
		\SetAlgoNoLine
		\SetKwBlock{enum}{typedef enum Status\{}{\}Status\;}
		\enum{OK=1,\\ ERROR=0}
	\end{algorithm}

	因为C语言中利用非$0$和$0$分别表示真与假，所以可以用枚举型变量中的两个值分别表示$0$和$1$。若程序运行过程中内存开辟失败，子函数立即停止运行并将$0$值返回，这样在主函数中我们就可以获知子函数是否能够正确运行，从而确定程序是否需要继续执行。为此，我们定义下面的枚举型变量来表示状态成功与否，之后的程序同样也用此结构表示动态内存空间开辟成功与否。同理程序运行过程中有时需要从文本文件中或二进制文件读取一些特定的参数，为确定程序是否成功访问到文件，我们同样用此结构表示文件的访问状态。
	
	\begin{algorithm}
		\SetAlgoNoLine
		\SetKwBlock{enum}{typedef enum Status\{}{\}Status\;}
		\enum{OK=1,\\ ERROR=0}
	\end{algorithm}
	
	子函数中包含如下算法，利用$malloc$分配内存失败返回$NULL$这一特性，我们很容易就可以获知内存分配是否成功。
	
	\begin{algorithm}
		\SetAlgoNoLine
		\SetKw{malloc}{malloc}
		\SetKw{pfloat}{(float*)}
		\SetKw{float}{float}
		\SetKw{size}{sizeof}
		\SetKwIF{mIf}{mElseIf}{mElse}{if}{}{else if}{else}{}
		ric= \pfloat \malloc ( \size ( \float )*nt )\;
		\mIf{(!ric)}{\Return{ERROR}\;}
		\Return{OK}\;
	\end{algorithm}
	
	最后，通过输入子波的主频和采样点的个数可以得到不同的子波图像。
	
	\begin{figure}[!h]
		\begin{tabular}{cc}
			\begin{minipage}{0.48\linewidth}
				\centerline{\includegraphics[width=6cm]{rick1.eps}}
			\end{minipage}			
			\begin{minipage}{0.48\linewidth}
				\centerline{\includegraphics[width=6cm]{rick2.eps}}
			\end{minipage}
		\end{tabular}
		\setlength{\belowcaptionskip}{2pt} 
		\caption{不同参数下的雷克子波图像}
	\end{figure}
	\newpage
	\section{傅里叶变换}
	\subsection{准备知识}
	傅里叶变换可以将满足一定条件的函数表示成正弦（余弦）函数的和。
	
	
	%%傅里叶变换
	\begin{equation}
		\label{fourier}
		\begin{aligned}
			& F\left(\omega\right) = \mathcal{F}\left[f\left(t\right)\right] = \int_{-\infty}^{\infty}f\left(t\right)e^{-i \omega t}\mathrm{dt}\\
			\\
			& f\left( t \right) = \mathcal{F}^{-1}\left[F\left(\omega\right)\right] =  \frac{1}{2\pi}\int_{-\infty}^{\infty}F\left(\omega\right)e^{i \omega t}\mathrm{d \omega}
		\end{aligned}
	\end{equation}

	为了使用计算机对信号进行傅里叶变换，通常对连续的信号进行等间距采样，采取有限数量的点，进行离散傅里叶变换。离散傅里叶变换的公式如下
	%%离散傅里叶变换
	
	\begin{equation}
		\label{fourier11}
		F\left(\omega\right)=\sum_{t=0}^{N-1}f\left(t\right)e^{-i\frac{2\pi \omega t}{N}}
	\end{equation}

	若采样率为$f_s$，$F[\omega]$表示频率为$\frac{\omega k_s}{N}$的相关度。
	
	我们也可以通过离散傅里叶逆变换将信号从频率域变换到时间域，离散傅里叶逆变换程序如下
	
	\begin{equation}
		\label{fourier12}
		f\left(t\right)=\sum_{t=0}^{N-1}F\left(\omega\right)e^{i\frac{2\pi \omega t}{N}}
	\end{equation}
	
	根据时间间隔$dt$，逆变换得到的$f(t)$表示时间为为$t \cdot dt$的信号强度。
	
	同理将表示信号强度原函数按两个维度分别等距划分，然后沿着一个维度做离散傅里叶变换后的结果作为新的函数，对新函数在沿另一个维度再做离散傅里叶变换，即可得到二维傅里叶变换的值。通过查询资料，我了解到二维傅里叶变换的结果可以用来描述图片的亮度。同样通过对两个维度分别进行一次离散傅里叶逆变换，可以得到二维离散傅里叶逆变换。二维离散傅里叶正逆变换公式如下
	
	\begin{equation}
		\label{fourier21}
		F\left(u,v\right)=\sum_{x=0}^{M-1}\sum_{y=0}^{N-1}f\left(x,y\right)e^{-i\frac{2\pi u x}{M}}e^{-i\frac{2\pi v y}{N}}
	\end{equation}
	
	\begin{equation}
		\label{fourier22}
		f\left(x,y\right)=\sum_{u=0}^{M-1}\sum_{v=0}^{N-1}F\left(u,v\right)e^{-i\frac{2\pi u x}{M}}e^{-i\frac{2\pi v y}{N}}
	\end{equation}
	
	\subsection{程序设计}
	通过查询相关知识，我们得到了式\ref{fourier11}到式\ref{fourier22}四个式子所表示的离散傅里叶变换公式。但即便有了公式，我们仍不能直接用C语言将该数学表达式描绘出来，其原因在于C语言中并没有复数这一数据结构。为了实现程序，我们需要定义复数这一数据类型及其相关算法。
	
	复数的数据类型如下，其中$r$和$i$分别表示复数的实部和虚部。
	\begin{algorithm}
		\SetAlgoNoLine
		\SetKwBlock{Complex}{typedef struct Complex\{}{\}Complex\;}
		\SetKw{float}{float}
		\Complex{\float r\;\float \i\;}
	\end{algorithm}
	
	在离散傅里叶正逆变换中，出现了复数之间的计算，包括复数的加法、乘法、幂运算等。为了提高代码的可读性，我将需要用到的复数运算法则封装成函数，供程序运行过程中调用。
	
	在傅里叶正逆变换中，因为加和符号的存在，我们都需要将得到的结果一一加起来，这就涉及到复数的加法运算。复数加法的算法定义如下：
	\begin{algorithm}
		\SetAlgoNoLine
		\SetKw{Complex}{Complex}
		\SetKwBlock{add}{Complex add(Complex z1,Complex z2)\{}{\}}
		\add{\Complex z\;z.r=z1.r+z2.r\;z.i=z1.i+z2.i\;\Return{z}\;}
	\end{algorithm}

	傅里叶逆变换中，$F$是一个复函数，积分的核是一个复数的指数，且这个复数的实部恒为0；傅里叶正变换中，$f$是一个实函数，积分的核是一个复数的指数,且这个复数的实部恒为0，因此需要分别定义复数和纯复指数乘积的运算法则以及实数和纯复指数乘积的运算法则。
	
	\begin{algorithm}
		\SetAlgoNoLine
		\SetKw{Complex}{Complex}
		\SetKw{float}{float}
		\SetKwBlock{multi}{Complex multi(Complex z1,Complex z2)\{}{\}}
		\multi{\Complex temp\;temp.r=z1.r*cos(z2.i)-z1.i*sin(z2.i)\;temp.i=z1.r*sin(z2.i)+z1.i*cos(z2.i)\;\Return{temp}\;}
	\end{algorithm}
	
	
	\begin{algorithm}
		\SetAlgoNoLine
		\SetKw{Complex}{Complex}
		\SetKw{float}{float}
		\SetKw{multi}{multi}
		\SetKw{real}{real}
		\SetKwBlock{real}{Complex real(float r1,float r2)\{}{\}}
		\real{\Complex z\;z.i=0\;z.r=r1\;z=\multi (z,r2)\;\Return{z}\;}
	\end{algorithm}
	
	\newpage
	
	
	根据式\ref{fourier11}以及定义的复数运算法则，我们定义如下算法表示一维离散傅里叶变换	
	\begin{algorithm}
		\SetAlgoNoLine
		\SetKw{pie}{define PI 3.1415926\\}
		\SetKwBlock{foro}{for(iw=0;iw<nw;iw++)\{}{\}}
		\SetKwBlock{fort}{for(it=0;it<nt;it++)\{}{\}}
		\SetKw{add}{add}
		\SetKw{PI}{PI}
		\SetKw{real}{real}
		\pie
		\foro{\fort{amp[iw]=\add (amp[iw],\real(rick[it],-2.* \PI *it*iw/nt))\;}}
	\end{algorithm} 

	其中$nt$表示一共有$nt$个采样点，$it$表示第$it$个采样点，对于时间采样间隔为$dt$的样本，其采样率为$\frac{1}{dt}$。$nw$表示变换后生成的点数，$iw$表示生成的第$iw$个点。
	
	根据式\ref{fourier12}以及定义的复数运算法则，我们定义如下算法表示一维离散傅里逆叶变换		
	\begin{algorithm}
		\SetAlgoNoLine
		\SetKw{pie}{define PI 3.1415926\\}
		\SetKwBlock{foro}{for(it=0;it<nt;it++)\{}{\}}
		\SetKwBlock{fort}{for(iw=0;iw<nw;iw++)\{}{\}}
		\SetKw{add}{add}
		\SetKw{PI}{PI}
		\SetKw{real}{real}
		\pie
		\foro{\fort{rick[it]=\add (rick[it],\real(amp[iw],2.* \PI *it*iw/nw)).r\;}}
	\end{algorithm}
	
	因为原信号为实信号，为减少算法的空间复杂度，此算法在计算过程中直接选取了逆变换后的实部。经验证，变换后的虚部在0附近徘徊，基本忽略不计。
	
	由以上分析和算法设计，我们可做出特定雷克子波的傅里叶变换图像，如下图所示。
			\begin{figure}[!h]
		\begin{center}
			\includegraphics[scale=0.2]{fourier.png}
			\caption{特定雷克子波的傅里叶变换}
		\end{center}
	\end{figure}
\newpage
	二维傅里叶变换的对象是播的反射过程，我们对波在时间和空间维度分别进行傅里叶变换。波的反射过程如下：在地面上某点放置一个波源，其向下面的各个方向连续发射雷克子波，子波在某一深度经反射传回地面，这样我们就可以在地面上的不同位置接收到子波。
	假设子波反射满足反射定律，我们可做出如下示意图。
	\begin{figure}[!h]
		\begin{center}
			\begin{tikzpicture}
				\draw (0,0)--(10,0);
				\draw (0,2)--(10,2);
				\fill (5,2) circle (2pt);
				\draw[-stealth] (5,2)--(4,0);
				\draw[-stealth] (4,0)--(3,2);
				\draw[-stealth] (5,2)--(3,0);
				\draw[-stealth] (3,0)--(1,2);
				\draw[-stealth] (5,2)--(6,0);
				\draw[-stealth] (6,0)--(7,2);
				\draw[-stealth] (5,2)--(7,0);
				\draw[-stealth] (7,0)--(9,2);
				\draw[|<->|] (-2,2)--(-2,0);
				\node (depth) at(-1,1) {depth};
				\node (sx) at(2.5,3.5) {sx};
				\draw[|<->|] (0,4)--(5,4);
			\end{tikzpicture}
		\end{center}
		\caption{波的反射示意图}
	\end{figure}
	
	根据式\ref{fourier21}和\ref{fourier22}再加上之前定义的复数运算法则，我们可以直接用循环结构将二维变换写出。但如果这样就要使用四层循环嵌套，时间复杂度为$\mathrm{O(M^2N^2)}$，对于时间长度上千的子波都要计算几个小时。如果先对每道子波沿时间维度进行一次傅里叶变换，在对得到的新函数沿空间维度进行一次傅里叶变换，则新算法的时间复杂度仅为$\mathrm{O(MN^2+NM^2)}$，尽管阶次仍然很高，但与之前相比，消耗的时间大大减少。
	
	经过如上分析，分别定义如下算法表示二维离散傅里叶正逆变换
	
	\begin{algorithm}
		\SetAlgoNoLine
		\SetKw{pie}{define PI 3.1415926\\}
		\SetKwBlock{foru}{for(iu=0;iu<nx;iu++)\{}{\}}
		\SetKwBlock{forv}{for(iv=0;iv<nv;iv++)\{}{\}}
		\SetKwBlock{fort}{for(it=0;it<nt;it++)\{}{\}}
		\SetKwBlock{forx}{for(ix=0;ix<nx;ix++)\{}{\}}
		\SetKw{add}{add}
		\SetKw{PI}{PI}
		\SetKw{real}{real}
		\SetKw{multi}{multi}
		\SetKw{DFT}{DFT}
		\pie
		\foru{\DFT(nt,dt,\&dv,rick[iu],\&amp[iu])\;}	
		\foru{\forv{amp1[iu][iv]=amp[iu][iv]\;}}
		\forv{\foru{\forx{amp[iu][iv]=\add(amp[iu][iv],\multi(amp1[ix][iv],-2*PI/nx*ix*iu))\;}}}
	\end{algorithm}

		\begin{algorithm}
		\SetAlgoNoLine
		\SetKw{pie}{define PI 3.1415926\\}
		\SetKwBlock{foru}{for(iu=0;iu<nx;iu++)\{}{\}}
		\SetKwBlock{forv}{for(iv=0;iv<nv;iv++)\{}{\}}
		\SetKwBlock{fort}{for(it=0;it<nt;it++)\{}{\}}
		\SetKwBlock{forx}{for(ix=0;ix<nx;ix++)\{}{\}}
		\SetKw{add}{add}
		\SetKw{IDFT}{IDFT}
		\SetKw{PI}{PI}
		\SetKw{real}{real}
		\SetKw{multi}{multi}
		\SetKw{DFT}{DFT}
		\pie
		\forx{\IDFT(nw,dw,\&dt,amp[ix],\&rick[ix])\;}	
		\forx{\fort{rick1[ix][it]=rick[ix][it]\;}}
		\fort{\forx{\foru{rick[ix][it]= \add (rick[iu][iv],\real(rick1[iu][iv],2*PI/nx*ix*iu)).r\;}}}
	\end{algorithm}
\newpage
经过如上分析和算法设计，我们做出反射波的二维傅里叶变换图像以及反变换的实部。

		\begin{figure}[!h]
			\begin{center}
		\includegraphics[scale=0.3]{newf.png}
		\caption{二维傅里叶变换}
	\end{center}
		\end{figure}

	\begin{figure}[!h]
		\begin{center}
			\includegraphics[scale=0.3]{data.png}
	\caption{二维逆变换的实部}
		\end{center}
	\end{figure}
	\section{波动方程数值模拟}
	\subsection{准备知识}
	波在介质中的传播与介质的密度和压强有关，通过查询相关知识，在波场中有如下方程成立。
	
	\begin{equation}
		\label{wave}
		\begin{aligned}
			\frac{\partial v_x}{\partial t} &= -\frac{1}{\rho}\frac{\partial p}{\partial x}\\
			\\
			\frac{\partial v_z}{\partial t} &= -\frac{1}{\rho}\frac{\partial p}{\partial z}\\
			\\
			\frac{\partial \rho}{\partial t} = -\rho & c^2 \left(\frac{\partial v_x}{\partial x}+\frac{\partial v_z}{\partial z} \right)
		\end{aligned}
	\end{equation}

	\subsection{程序设计}
	波场模拟依据式\ref{wave}，但在C语言中无法进行微分和积分运算，我们可以用差分和求和分别代表微分积分运算。

	\begin{equation}
		\label{diff}
		\begin{aligned}
			\frac{v_x(t+\Delta t)-v_x(t)}{\Delta t} &= -\frac{1}{\rho}\frac{p(x+\Delta x)-p(x)}{\Delta x}\\
			\\
			\frac{v_z(t+\Delta t)-v_z(t)}{\Delta t} &= -\frac{1}{\rho}\frac{p(z+\Delta z)-p(z)}{\Delta z}\\
			\\
			\frac{p\left(t+\Delta t\right)-p\left(t\right)}{\Delta t}=-\rho c^2 &\left(\frac{v_x\left(x+ \Delta x\right)-v_x\left(x\right)}{\Delta x}+\frac{v_z\left(z+ \Delta z\right)-v_z\left(z\right)}{\Delta z}\right)
		\end{aligned}
	\end{equation}
	
	根据以上分析，定义应力算法如下\newpage
	\begin{algorithm}
		\SetAlgoNoLine
		\SetKwBlock{forx}{for(ix=1;iu<nx-1;iu++)\{}{\}}
		\SetKwBlock{forz}{for(iz=1;iz<nz-1;iz++)\{}{\}}
		dx0 = dt / dx\;
		dz0 = dt / dz\;
		\forx{\forz{p[ix][iz] += ((vx[ix][iz] - vx[ix - 1][iz]) * dx0 + (vz[ix][iz] - vz[ix][iz - 1]) * dz0) * den[ix][iz] *v[ix][iz] * v[ix][iz] + data[ix][iz]\;}}
	\end{algorithm}
	
	波速对应的算法如下
		\begin{algorithm}
		\SetAlgoNoLine
		\SetKwBlock{forx}{for(ix=0;iu<nx-1;iu++)\{}{\}}
		\SetKwBlock{forz}{for(iz=0;iz<nz-1;iz++)\{}{\}}
		dx0 = dt / dx\;
		dz0 = dt / dz\;
		\forx{\forz{vx[ix][iz] += 1. / den[ix][iz] * (p[ix + 1][iz] - p[ix][iz]) * dx0\;vz[ix][iz] += 1. / den[ix][iz] * (p[ix][iz + 1] - p[ix][iz]) * dz0\;}}
		\end{algorithm}

	因为波速和应力相互影响，因此可以采用网格法分别对应力和波速依次求解,依次调用应力算法和波速算法即可求解出最后时刻对于每一个坐标点的波场参数。
	
	为了更加贴近物探的实际应用，本程序的介质取为如下图所示的正断层。
	\begin{center}
		\begin{tikzpicture}
			\draw[scale=0.5] (0,0)--(-4,4)--(-8,4)--(-8,0)--cycle;
			\draw[scale=0.5] (-2,2)--(2,-2)--(6,-2)--(6,2)--cycle;
		\end{tikzpicture}
	\end{center}

在断层中起始波速为$1500m/s$，为了简化模型，只分析断层下盘底部到顶部这一高度的波场,假s上盘以上部分的波速取为$700m/s$，以此建立波速矩阵，代表初始波速。波速矩阵如下所示

	$$
	V=
	\begin{bmatrix}
		1500 & 700  & \cdots   & 700   \\
		1500 & 1500  & \cdots   & 700  \\
		\vdots & \vdots  & \ddots   & \vdots  \\
		1500 & 1500  & \cdots\  & 1500  \\
	\end{bmatrix}
	$$
	
	交错调用应力和速度算法，即可模拟出该正断层中的波场。最后时刻波场中$p$的分布如下图所示
	\begin{figure}[!h]
		\begin{center}
			\includegraphics[scale=0.3]{wave.png}
			\caption{断层波场模拟}
		\end{center}
	\end{figure}
	
	\section{实训总结及感悟}
	此次实训，我不仅重新复习了C语言的知识点，又学习到了动态分配的知识。老师在课堂上讲解的十分详细，尽管实训过程中遇到了许多没有学过的知识点，老师也耐心给我们讲解，让我们学到了一些地震波的相关知识。
	
	课上老师会先带着我们浏览一下程序和目标图像，我在课上先跟着老师写，把每一个参数的含义想清楚弄明白，最后再查找相关资料，整合后写出自己的程序。例如在离散傅里叶变换这一节，我没有在程序之初便指定频率间隔这一参数，而是利用采样率和采样点的个数来计算出频率间隔。
	
	因为我的电脑操作系统为Linux，不能使用Fimage对数据进行分析，因此在画图的时候我抽出一些时间学习了Seismic Unix，我相信对之后物探地震分析也有很大的帮助。
	
	通过实训，我独立编写程序的能力得到了提高，能够较熟练地使用计算机编程解决实际问题。最后感谢学习、感谢老师给我们这次动手的机会，让我们可以亲手利用计算机解决实际工程中的一些问题。
	\section{实训日志}
	\month=6
	\noindent
	\number\month .\thedirary
	
	今天是第一天实训，今天实训的任务是人为输入参数，利用C语言计算出给定参数的雷克子波信号数据，并利用EXCEL将其信号图画出。课上，老师首先给我们介绍了雷克子波及其模型,然后为我们演示了程序,最后我们自己写程序,通过设置参数得到了不同的子波数据,并利用EXCEL绘制出波形.
	
	\addtocounter{dirary}{1}
	\noindent
	\number\month .\thedirary
	
	今天是实训的第二天,今天的任务是编写程序完成离散傅里叶正变换和逆变换,并对前一天计算出的雷克子波进行傅里叶变换，对变换得到的数据再进行反变换，比较前后数据是否相同。
	
	\addtocounter{dirary}{1}
	\noindent
	\number\month .\thedirary
	
	今天是实训的第三天，今天的主要任务是将之前写的傅里叶变换封装成子函数，并设计程序计算二维傅里叶变换。今天的程序以雷克子波的反射为例，进行二维傅里叶正反变换。
	
	\addtocounter{dirary}{1}
	\noindent
	\number\month .\thedirary
	
	今天是实训的第四天，今天的主要任务是将二维傅里叶正反变换设计成子函数，并验证是否正确。
	
	\addtocounter{dirary}{1}
	\noindent
	\number\month .\thedirary
	
	今天是实训的第五天，今天主要是完成昨天未完成的程序，继续将二维傅里叶变换编写完整。
	
	\addtocounter{dirary}{1}
	\noindent
	\number\month .\thedirary
	
	今天是实训的第六天，今天主要是自我编程练习。我将之前写的一维傅里叶变换函数进行了改写，增加了判断动态分配是否成功的返回值，并将复数的一些相关运算封装成函数，方便复用。
	
	\addtocounter{dirary}{1}
	\noindent
	\number\month .\thedirary
	
	今天是实训的第七天，主要是自己编程练习。之前写的二维傅里叶变换仍未完成，我查找了相关资料，根据老师上课讲解的例子继续完成程序。
	
	\addtocounter{dirary}{1}
	\noindent
	\number\month .\thedirary
	
	今天是实训的第八天，今天的任务是完成波动方程数值模拟的编程。通过之前数理方程的学习，我对波动方程有了一定的认识，老师引导我们使用差分法来进行程序设计。
	
	\addtocounter{dirary}{1}
	\noindent
	\number\month .\thedirary
	
	今天的主要任务是继续完成前八天未完成的程序设计。我完成了昨天为写完的波动方程模拟程序，并绘制出来波场图像。
	
	\addtocounter{dirary}{1}
	\noindent
	\number\month .\thedirary
	
	今天是实训的最后一天，我先是对之前未完成的二维傅里叶变换进行了完善。之后又对波场相关的参数进行了设置修改，以贴合断层中波的传递。
	
	\newpage
	\begin{appendices}
		\dottedcontents{section}[4em]{\normalsize}{4.0em}{4pt}  %%避免文字重叠
		\titleformat{\section}{}{\thesection}{2em}{}
		\renewcommand{\thesection}{附录\Alph{section}.}
		\section{雷克子波的生成}
		\begin{lstlisting}
Status wave(float fpeak,float dt,int n,float **amp)
{
	float temp;
	float T = 1. / fpeak;
	float t;
	int it;
	*amp = (float *)malloc(sizeof(float) * n);
	if(!(*amp))
		return ERROR;
	for (it = 0; it < n;it++)
	{
		t = it * dt - T;
		temp = PI * fpeak * t;
		temp = temp * temp;
		(*amp)[it] = (1 - 2 * temp) * exp(-temp);
	}
	return OK;
}
		\end{lstlisting}
		\section{一维傅里叶变换}
		\begin{lstlisting}
Status DFT(int nt,float dt,float *dw,float *rick,Complex **amp)
{
	int iw,it;
	*dw=1./nt/dt;
	*amp=(Complex *)malloc(sizeof(Complex)*nt);
	for(iw=0;iw<nt;it++) (*amp)[iw].r=(*amp)[iw].i=0;
	for(iw=0;iw<nt;iw++)
	for(it=0;it<nt;it++)
	{
	(*amp)[iw]=add((*amp)[iw],real(rick[it],-2.*PI*it*iw/nt));
	}
	return OK;
}

Status IDFT(int nw,float dw,float *dt,float **rick,Complex *amp)
{
	int iw,it;
	*rick=(float *)malloc(sizeof(float)*nt);
	for(it=0;it<nt;it++) (*rick)[it]=0;
	for(it=0;it<nt;it++)
	for(iw=0;iw<nw;iw++)
	{
	(*rick)[it]+=1./nw*multi(amp[iw],2.*PI/nw*iw*it).r;
	}
	*dt=1./nw/dw;
	return OK;
}
		\end{lstlisting}
		\section{二维傅里叶变换}
		\begin{lstlisting}
//假设出现的所有变量均初始化为0，nx、nv等变量从文件中读取，u、v分别
//代表变换后的两个自变量

//正变换核心代码
for(iu=0;iu<nx;iu++){
	DFT (nt,dt,&dv,rick[iu],&amp[iu]);
}
for(iu=0;iu<nx;iu++){
	for(iv=0;iv<nv;iv++){
		amp1[iu][iv]=amp[iu][iv];
	}
}
for(iv=0;iv<nv;iv++){
	for(iu=0;iu<nx;iu++){
		for(ix=0;ix<nx;ix++){
			amp[iu][iv]=add (amp[iu][iv],multi (amp1[ix][iv],
			-2*PI/nx*ix*iu));
		}
	}
}

//逆变换核心代码
define PI 3.1415926
for(ix=0;ix<nx;ix++){
	IDFT (nw,dw,&dt,amp[ix],&rick[ix]);
}
for(ix=0;ix<nx;ix++){
	for(it=0;it<nt;it++){
		rick1[ix][it]=rick[ix][it];
	}
}
for(it=0;it<nt;it++){
	for(ix=0;ix<nx;ix++){
		for(iu=0;iu<nx;iu++){
			rick[ix][it]= add (rick[iu][iv],real (rick1[iu][iv],
			2*PI/nx*ix*iu)).r;
		}
	}
}
		\end{lstlisting}
		\section{波动方程数值模拟}
		\begin{lstlisting}
void updatep(float **p, float **vx, float **vz, float **v, float **den, float dx, float dz, float dt, int nx, int nz,
float **data) {
	int ix, iz;
	float dx0, dz0;
	dx0 = dt / dx;
	dz0 = dt / dz;
	for (ix = 1; ix < nx - 1; ix++)
	for (iz = 1; iz < nz - 1; iz++) {
		p[ix][iz] += ((vx[ix][iz] - vx[ix - 1][iz]) * dx0 +
		 (vz[ix][iz] - vz[ix][iz - 1]) * dz0) * 
		 den[ix][iz] *v[ix][iz] * v[ix][iz] + data[ix][iz];
	}
}

void updatev(float **p, float **vx, float **vz, float **den, float dx, float dz, float dt, int nx, int nz) {
	int ix, iz;
	float dx0, dz0;
	dx0 = dt / dx;
	dz0 = dt / dz;
	for (ix = 1; ix < nx - 1; ix++)
	for (iz = 1; iz < nz - 1; iz++) {
		vx[ix][iz] += 1. / den[ix][iz] * 
		(p[ix + 1][iz] - p[ix][iz]) * dx0;
		vz[ix][iz] += 1. / den[ix][iz] * 
		(p[ix][iz + 1] - p[ix][iz]) * dz0;
	}
}

int main() {
	FILE *fpin, *fpout;
	char out[100], in[100];
	int sx0, sz0, sd0;
	float fpeak, dt, *ric, dz, dx, **den;
	float sx, sdepth, gdepth;
	float **data, **p, **vx, **vz, **v, **record;
	int it, nt, ix, iz, nx, nz;
	fpin = fopen("../par/wave.txt", "r");
	fscanf(fpin, "fpeak=%f\n", &fpeak);
	fscanf(fpin, "dt=%f\n", &dt);
	fscanf(fpin, "nt=%d\n", &nt);
	fscanf(fpin, "out=%s\n", out);
	fscanf(fpin, "nx=%d\n", &nx);
	fscanf(fpin, "nz=%d\n", &nz);
	fscanf(fpin, "dx=%f\n", &dx);
	fscanf(fpin, "dz=%f\n", &dz);
	fscanf(fpin, "sx=%f\n", &sx);
	fscanf(fpin, "sdepth=%f\n", &sdepth);
	fscanf(fpin, "gdepth=%f\n", &gdepth);
	fclose(fpin);
	ric = rick(fpeak, dt, nt);
	sx0 = (int) (sx / dx + 0.5);
	sz0 = (int) (sdepth / dz + 0.5);
	sd0 = (int) (gdepth / dz + 0.3);
	p = (float **) malloc(sizeof(float *) * nx);
	vx = (float **) malloc(sizeof(float *) * nx);
	vz = (float **) malloc(sizeof(float *) * nx);
	v = (float **) malloc(sizeof(float *) * nx);
	den = (float **) malloc(sizeof(float *) * nx);
	data = (float **) malloc(sizeof(float *) * nx);
	record = (float **) malloc(sizeof(float *) * nx);
	for (ix = 0; ix < nx; ix++) {
		record[ix] = (float *) malloc(sizeof(float) * nt);
		p[ix] = (float *) malloc(sizeof(float) * nz);
		vx[ix] = (float *) malloc(sizeof(float) * nz);
		vz[ix] = (float *) malloc(sizeof(float) * nz);
		data[ix] = (float *) malloc(sizeof(float) * nz);
		v[ix] = (float *) malloc(sizeof(float) * nz);
		den[ix] = (float *) malloc(sizeof(float) * nz);
		for (iz = 0; iz < nz; iz++) {
			p[ix][iz] = 0;
			vx[ix][iz] = 0;
			vz[ix][iz] = 0;
			data[ix][iz] = 0;
			den[ix][iz] = 1;
			v[ix][iz] = 1500;
		}
	}
	
	for(iz=0;iz<nz;iz++)
	{
		for(ix=0;ix<nx;ix++)
		{
			if(iz>500)
			v[ix][iz]=1500;
			else if(ix-iz>=500)
			v[ix][iz]=700;
		}
	}
	
	fpout = fopen("../data/wave_rick.txt", "w");
	for (it = 0; it < nt; it++)
	fprintf(fpout, "%f\t%f\n", it * dt, ric[it]);
	fclose(fpout);
	for (it = 0; it < nt; it++) {
		data[sx0][sd0] = ric[it];
		updatep(p, vx, vz, v, den, dx, dz, dt, nx, nz, data);
		updatev(p, vx, vz, den, dx, dz, dt, nx, nz);
		for (ix = 0; ix < nx; ix++)
		record[ix][it] = p[ix][sd0];
		if (it % 4 == 0)
		printf("%d\n", it);
	}
	fpout = fopen("../data/wave.bin", "wb");
	for (ix = 0; ix < nx; ix++)
	fwrite(p[ix], sizeof(float), nz, fpout);
	fclose(fpout);
	fpout = fopen("../data/record.bin", "wb");
	for (ix = 0; ix < nx; ix++)
	fwrite(p[ix], sizeof(float), nt, fpout);
	fclose(fpout);
	for (ix = 0; ix < nx; ix++) {
		free(v[ix]);
		free(vx[ix]);
		free(vz[ix]);
		free(data[ix]);
		free(den[ix]);
	}
	free(data);
	free(den);
	free(vx);
	free(vz);
	free(v);
	free(ric);
	return 0;
}

		\end{lstlisting}
	\end{appendices}
\end{document}